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Abstract. We study one-dimensional conservation law. We use a character- 
istic surface to define a class of functions, within which the integral version of 
the conservation law is solved in a simple and direct way, avoiding the use of 
the concept of weak solutions. We develop a simple numerical method for com- 
puting the unique solution. The method uses the equal-area principle and gives 
the solution for any given time directly. 
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I. Introduction 

We study the one-dimensional scalar conservation law 

(1) u t + f(u) x = 0, 

where u : R x [0, oo) — > R represents the conserved quantity (such as mass, density, 
momentum) while / : R — > R is the flux. We equip the equation with the 
initial condition 

(2) u{x,0) = h{x), xeR. 

It is well known that classical solutions of Q-([2]), even for smooth initial condi- 
tions, do not always exist. Usually, the concept of solutions is generalized to the 
so-called weak solutions which are not unique and, in order to obtain the proper 
solution, one has to impose some entropy condition. There is a rich mathematical 
theory on this topic, see for example monographs [U [5j (6J [7J |8] . 

We proceed more directly. In order to allow discontinuities we consider the integral 
form of the conservation law 



d 



i, 



(3) j t J u(x, t)dx = f (u(a, £)) - f (u(b, t)) 

for all a < b and t > for which u is continuous in points (a, t) and (b, t). It says 
that the area J a u(x, t) dx changes in time according to the flux at the boundaries. 
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If both u and / are continuously differentiable, ^ implies ([TJ. We search for the 
solutions of ^ inside a class of functions T that is defined by a characteristic 
surface associated to the problem ([T])-([2]). This functions may have jumps along 



some locally smooth curves in the (x,t)-plane (see Definition 2.2). In Section[2]we 
prove uniqueness of solutions to ^ within the class T. 

Taking the integral over all R, Q implies 

(4) / u(x, t) dx= / h(x) dx for all t > 0. 

Jr Jr 

Hence the area under the graph of the solution does not change in time. Based 
on this observation we propose in Section [3] a simple method which we call an 
equal-area method. We show that, under suitable conditions, it yields the unique 
solution of (|3| within the class T. 

The idea of using the equal-area principle is not new. It has already been exploited 
in the classical textbook by Whitham [121 Sect. 2.8-2.9]. However, the method 
there is developed only analytically in a complicated way that is not suitable for 
explicit computations. Recently, Farjoun and Seibold |3j suggested a conservative 
particle method that uses equal-area principle. They use the Lagrangean approach, 
representing the solution as a cloud of particles which move with the flow. Particles 
carry function values and move according to their characteristic velocities. When 
the characteristic curves collide, the particles are merged in such a way that the 
total area under the function is conserved. So far we are not aware of any other 
numeric method using thr equal-area principle. 

In Section [4] we describe a numerical algorithm we used to implement the method. 
We discuss numerical results obtained by our method in Section [5] and compare 
it to the finite volume method Clawpack [2] as well as to the newly suggested 
conservative particles method Particleclaw [9]. 

Our method performs well. It is exact and the shocks are formed when necessary, 
with high accuracy. In contrast to other known methods, where the solution at 
selected time is obtained by evolution from initial conditions, it gives the solution 
for any given time directly. 



2. Characteristic surface and integral solutions 

A well-known method for treating the initial value problem ([l])-([2]) is the method 
of characteristics. The characteristics are curves in the (x, t)-plane along which 
the function u is constant. In our case they are lines of the form 



(5) 
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see e.g. (HI Sec. 2.2]. It is easy to see that, if a C 1 -solution u(x,t) of ([T])-([2]) exists, 
its graph in IR 3 is given by 



(6) 



r := {(x, t, y) = (£ + /' (HO) t, t, MO) |(ei,i>o}. 



We shall call T the characteristic surface to the problem (JXJ)- ([2]) . An example of 



such a surface is seen in Figure 1(a) 





(a) Parallel characteristics. (b) Characteristics collide. 

Figure 1. Examples of characteristic surfaces T. 

We can form T a-priori, before investigating the solvability of ([lJ-Q, whereby 
it might happen that it represents a multivalued function which cannot be the 
solution of our problem (see Figure 1(b)). Indeed, this problem occurs whenever 



the characteristics collide in the (x, t)-plane. In this case the proper solution has 
jumps along some curves in the (x, i)-plane. Its graph, however, is still a subset of 
T. Starting with T, the solution can thus be obtained by finding the appropriate 
position of the jumps. 

In the case when the initial function h is not continuous the above defined char- 
acteristic surface T is not connected. Before proceeding we shall hence modify the 
definition of characteristic surface to correct this. For some fixed t > we define 
plane transformation Gt as 



(7) 



Gt(x,y) := (x + f'(y)t,y). 



Denote by 70 the graph of the initial function h together with vertical line segments 
joining discontinuities and 



7t := G t (7o) 



which is continuous curve for all t > 0, see Figure [2] 
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Figure 2. The graph of h and the curves 70 and 7 f . 



Definition 2.1. The bounding characteristic surface to the problem ([T])-([2]) is 
defined as 



Note that T contains the characteristic surface T and that the two surfaces agree 
whenever h is continuous. 

We are now ready to define the appropriate domicile for our solutions. 

Definition 2.2. Let / G C 2 (1R) and h is a piecewise C 1 -function with compact 
support. We say that u G T = T(/, h) if the following conditions are satisfied: 

(i) the function u = u(x,t) is defined everywhere on R x [0, 00), 

(ii) the function u(x, 0) = h(x) for all xel, 

(iii) the graph of u in IR 3 is a subset of the bounding characteristic surface Y 
defined by / and h, 

(iv) the boundary of the graph of a is a finite union of C 1 -curves, and the 
projections of these curves on (x, t)-plane intersect any line t = t only 
finitely many times, 

(v) the integral u(x,t) alx is a continuous function of t for t > 0. 

Our aim is to find the solutions to the integral form of the conservation law (J3]) 
within the class T. First we demonstrate that every u G T solves Q on its 
of smoothness. 

Lemma 2.3. Let u G T be a C l -function on an open set D C R x [0, 00). Then 
u is a solution of Q on D. 

Proof. If the graph of u(D) lies on the original characteristic surface T, then u 
solves Q on D by the method of characteristics (see e.g. [5], §3]). Therefore we 
may assume that the graph of u(D) is contained in T \ T. 

We will show that the added vertical lines and their convolutions in time also yield 
a solution in a similar way as the original characteristics do. For any fixed £ we 



(9) 



r := {(%, t, y) I (x, y) G 7*, £ > 0}. 
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parametrize the complemented surface in f \ T as 

(x,t,u(x,t)) = ^ + f'(r)t,t,T). 

Now by implicit derivation of the x-coordinate, by the equality r = r(x,t) = 
u(x,t), and by the Implicit function theorem (see [101 Theorem 9.28]) one obtains 

„,= » and „- ''< T > 



l"{r)t l"{r)f 
and it is easy to see that w, implicitly defined this way, is a local solution of ([!]). □ 

If, however, «GTis not smooth for some t > and i6R, then it has discontinu- 
ities called shocks which are positioned along piecewise smooth curves x = s(t) in 
the (x, t)-plane, called shock paths. By 2.2 iv), there are finitely many shock paths, 
which are all locally smooth. These paths may have singular points, they may cross 
or collide, but we can exclude these singularities without loss of generality. 

We now continue our treatment by showing that the well-known Rankine-Hugoniot 
condition holds in our setting. 

Lemma 2.4. A function u £ T is a solution of (|3| if and only if the following 
Rankine-Hugoniot condition is satisfied at all the shocks: 

n , f(u + ) - f(vr) 
(10) s'(t) = 



u + — u 



(by u + and u we denote the one-sided limits of the solution u(x,t) from the left 
and from the right side of the shock, respectively, i. e. u + (x ,t) = ]xni x \^ Bo u(x,t)). 

Proof. For the sake of simplicity we shall omit the arguments of functions whenever 
they are clear from the context. Take any a < b. If u 6 T is smooth for x G (a, b) 
and t > 0, it solves Q on (a, b) and we have 

d f b f b f b 

u dx = u t dx = — / f(u) x dx = f {u(a)) — f {u(b)) . 



dt 

Now assume that u has a shock on (a, b) x {t} with shock path x = s(t). From 



the condition (iv) of Definition 2.2 it follows that x = s(t) is a locally C 1 -function. 
Hence, we can compute 

d [ b J d [ a ® J d [ b J 
— u dx = — I u dx + — / u dx 
dt J a dt J a dt J s{t) 

rs(t) /■& 

(11) = / Ut dx + s'{i)u~ + j ut dx — s'(t)u + 

Ja Js(t) 

= f (u(a)) - f («") + / («+) - / («(&)) + s'(t) (u- - n + ) . 

Thus, u is a solution of ^ on (a, b) if and only if the Rankine-Hugoniot condition 
( 10 ) is fulfilled at the shock. If u has more then one shock on (a, b) x {t}, we divide 



(3 



MARJETA KRAMAR FIJAVZ, MITJA LAKNER, MARJETA SKAPIN RUGELJ 



the interval according to the shocks and proceed in the same manner. Note that 
Definition 2.2[ iv) implies that u has only finitely many shocks. □ 

We can prove uniqueness of the solutions to ([2])-((3]) in the class T, provided that 
the flux is a concave function and the solutions only have jumps upwards. The 
same is true for convex flux and downwards jumps. 

Proposition 2.5. Let f G C 2 (1R) be strictly concave. Letu,v G T be two solutions 
of (J3]) with the properties 

(a) u~ < u + and v~ < v + at every shock for u and v, respectively, and 

(b ) u(x, 0) = v(x, 0) for all x G R. 

Then the L l -norm 

(12) IK',*) ~v(; = forallt>0. 



Proof. We slightly modify the idea Lax used for proving [51 Theorem 3.4]. Let 
u, v G T be two solutions to ^ with u(x, 0) = v (x, 0) = h(x) for all x G K. Since 
they both lie on the same characteristic surface their difference u — v can change 
the sign only in the shocks of either (or both) of these two functions. Denote these 
sign-changing curves in (x,t) -plane by yk{t) and order them as 

yi(t) <V2(t) < ■■■ < y n +i(t). 

We work on maximal open intervals of t where the curves yk{t) do not intersect 
and the number of these curves does not change. For almost every t > we can 
then write 

rVk+i(t) 

(13) \\u(-,t) - v(-,t)\\i = y^pfc(t) where p k {t) = / \u(x,t) - v(x,t)\ dx. 

k=i Jyhit) 



Now choose any interval (yk(t),yk+i(t)), 1 < k < n. Without loss of generality we 
may assume that u(x,t) > v(x,t) on this interval, hence the absolute value under 
the integral defining pk{t) can be omitted. 

Note that Pk{t) is a (continuous) piecewise different iable function of t. If either u 
or v has some shocks inside the interval, in each of them the Rankine-Hugoniot 
condition is satisfied by Lemma 2A Dividing the interval according to these shocks 
and computing the derivative according to this division, similarly as it was done 
in (11), we see that the values around the shocks cancel out and only the values 
in yk{t) and yk+i{t) are important. Therefore we shall from now on assume that 
none of u and v has shocks inside the interval (yk(t),y k+1 (t)). 

Thus we may assume that function pk{t) is differentiable and we will now compute 
its derivative. As usual we will omit the arguments of functions whenever possible. 
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Since u and v solve ([TJ inside the interval and since the shock paths are piecewise 
differentiable we have 
rVk+i 

P'k(t) = ( u t~ v t ) dx + (u~(y k+ i) - v~{y k+1 )) y' k+1 - (u + (y k ) - v + {y k )) y' k 

■JVk 

(14) = / {u + (y k )) - / (v + (y k )) - (u + (y k ) - v + (y k )) y' k 

(15) - [/ (vT(y k +i)) - f (v~{y k +i)) - (vr{y k+1 ) - u~(y fc+1 )) y' k+1 ] ■ 

From now on we shall explain the calculations only for the left endpoint of the 
interval, the right endpoint can be treated symmetrically. By assumption, u > v 
inside the interval, hence u — v changes sign at the endpoints and all the jumps 
are upwards. Taking all this into account we see that u has a shock in y k while 
v may have a shock or not (in the latter case we take v~ = v + ). Furthermore we 
have 

(16) u~ < v~ < v + < u + . 

Applying the Rankine-Hugoniot condition for the speed of shock y' k for u we see 
that (14) equals 

7(« + )-/(« + ) f(u + )-f(u-y 



u 



+ _ 



Vu — u~ 



{u + - v + ) < 0, 



since by concavity of / and condition ( 16 ) the first factor is smaller or equal to 
while u + — v + > 0. 

Following the same line of arguments for the right endpoint y k +i we obtain the 
same conclusion for (15). 



We have thus shown that p' k (t) < 0, 1 < k < n. By (13), —v(- : t)\\i is then 

a decreasing function of t. Since by assumption ||it(-,0) — u(-,0)||i = we finally 
obtain (p}. □ 



Note that the condition (a) of Proposition 2.5 is actually an entropy condition (see 
e. g. [6, p. 36]). 



3. The equal-area method 



We now describe the equal-area method for obtaining the solutions starting from 
the bounding characteristic surface T defined in ([9]). First note that the trans- 
formation G t defined in ^ preserves area, since its Jacobian equals 1. Hence 
the area bounded by the curves jt and the x-axis remains unchanged in time and 
equals the initial area given by h(x) dx (compare the shaded areas in Figure 

\2\j. Intersecting f with the plane t = t for some fixed time to yields 7 <0 defined in 
^8]), see also Figure [2j Our strategy is to insert vertical cuts to 7j in such a way 
that the areas of the cut-off lobes coincide. Thus the initial area will be preserved. 



8 



MARJETA KRAMAR FIJAVZ, MITJA LAKNER, MARJETA SKAPIN RUGELJ 



Carrying out this procedure for all t we obtain a bounded, piecewise continuous 
function u(x, tY, whose graph, without the added vertical surfaces, is contained in 
f, see Figures and [i] 





(a) The bounding characteristic surface T. 



(b) The solution u(x, t) (with vertical sur- 
faces at the positions of the jumps). 



Figure 3. Graphs for the problem (JXJ)— (|2]) for the initial function 
given in (25). 



Lemma 3.1. Let f e C 2 (IR) and let h be a piecewise C l -function with compact 
support. Then for the solution u obtained by the above described equal-area method 
the following holds. 



(a) At every shock u satisfies the Rankine-Hugoniot condition (10) and has 
only jumps upwards: u~ < u + . 

(b) The shock paths are piecewise C 1 -curves. 

Proof. For t > we have the curves 

For £1 < £2 we close the "S'-curve" 7t([£i, £2]) with line segment between endpoints 
7t(£i) and 7* (£2); see Figure |lj The signed area defined by this closed curve by 
Green's Theorem equals 



1 ^ 



(17) + - / u x t (eo + a - (wtceo - y t m 







- yt (eo + (1 - j/t(£ 2 )) (x t (6) - ^(6))] 
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!) 




We define the mapping F : M 2 x IR+ — y M 2 by 

F(£i,&,t) :=(pt(£i,&W£i) -**(&)). 

In points where 

(18) xt{&)=*t{&), 

the determinant of the 2x2 Jacobian matrix ^ is equal to 

and is nonzero in jumps given by our method. 



If in addition to (18), we have the equal-area condition Pt„(£i, £2) = 0, then the 
Implicit function theorem (see [10, Theorem 9.28]) gives the existence of two C l - 
functions £2(^)5 such that for all t in some neighborhood of t Q the following 

holds: 

F(Ut)M),t) = (0,0), Uto) = £?, §(*o) = £ 2 °- 

This means that the equal-area condition holds for all parameters t in this neigh- 
borhood. Moreover, the Implicit function theorem gives us the formula 



ei(*o) 
8 (to) 



dF 



«9F 



Using symbolic computation (Mathematica [13]) one obtains 

/(*«?)) - f(He 2 )) + M) - 



£i(*c 



(^)-^))(l + / ff W))^)*o) 



By above, the shock path 



*(*) = &(*) +/'(%(*)))* 
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is locally a C 1 -function, which proves (b). Differentiating this function at the point 



t = to, we finally get the Rankine-Hugoniot condition (10) 

x'(t ) = g(to) + f"(h(£))ti(g)l(to)t + f(h($)) 



(19) 



u — u~ 



In the case when / is concave, /' is a decreasing function and the curves j t are for 
t > inclined to the left (regarding x-axis). It then follows from the construction 
that the solution u only has jumps upwards, hence also the assertion (a) is proved. 

□ 

We have proved the local smoothness of the shock paths, a nice property which is 
often assumed in advance and rarely verified. Assuming a finite number of shocks 
we are finally able to prove that the equal-area method yields the unique solution 
of our problem. 

Theorem 3.2. Let f e C 2 (R) be a strictly concave function and h a piecewise 
C 1 -function with compact support such that the function 

(20) x t {Z) = Z + f'(h(Z))t, £eM 

has only finitely many local extrema for every t > 0. Then the above described 
equal-area method yields the unique solution u £ T to S. 



Proof. First observe that conditions (i), (ii), (iii), and (v) of Definition 2.2 are 
trivially fulfilled for a function u obtained by the equal-area method. Since the 
number of local extrema of the function (|20l) is finite, so is the number of shocks 



obtained by the equal-area method. Hence the condition (iv) of Definition 2.2 



is 



satisfied by Lemma 3.1. Moreover, using Lemmas 3.1, 2.3, 2.4, and 2.5 we see that 
u E T is the unique solution of (j3]). □ 

A brief comment is in order here. For our theoretical approach we need to as- 
sume that there are only finitely many shocks. We are not aware of any explicit 
conditions in terms of functions / and h to meet this assumption (some generic 
conditions are given in [U HI]). In praxis however, it is not difficult to check the 



finiteness of the number of local extrema of the function Xt(£) in (20) for any 
given /, h, and t (note that our procedure gives a solution for any fixed time £!). 
Moreover, numerically this condition is always satisfied, since we use polygonal 
approximation for continuous curves. 
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4. The algorithm 

We shall now describe the algorithm based on the procedure introduced in the 
previous section. The solution at some fixed time to is obtained directly, without 
need to proceed in time. 

We start by taking a polygonal approximation Kq for the continuous curve 
(21) 7 t :=^ (7o), 

see ([7])-@. Traveling along the curve we gradually 'equalize' the areas. The 
obtained graph of solution is a subset of K (see Figure [6]). 

Iterative Step of the Algorithm 4.1. Let the parametrization (x(r), y(r)), r£l, 
of the curve Ki on the i-th step be such that x(t) is increasing on the far ends 
of the interval (— oo, oo). First we define three significant points for the curve Ki 
(see Figure [5]): 

(1) (5 = x(ti) is the first local maximum of x(t) in the direction of the in- 
creasing parameter r. If such a maximum does not exist, we are done and 
Ki is the graph of the weak solution (with redundant vertical lines in the 
jumps). 

(2) a = x(t 2 ) is is the first local minimum of x(t) from T\ onwards. Since the 
function x(t) in not bounded from above, such a minimum always exists. 

(3) 7 is the minimum of /3 and the first next local maximum of x(t). If such a 
maximum does not exist, let 7 = (3. 




Figure 5. Significant points a, f3, 7, and 5 on the curve Ki and 
the areas P\{$) and ^2(7)- 



Now we compute the areas. For any x G (a, (3) denote by pi(xo) the area bounded 
by the line x = xq and the part of the curve Ki that contains (x(ti), y(ri)). For 
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any xq G (a, 7) denote by ^2(^0) the area bounded by the line x 
of the curve Ki that contains (x(r 2 ), £/(t 2 )). For x G (a, 7) let 



Xq and the part 



(22) 



p(x) 



p 2 (x) -pi(x) 



Then pi(x) is continuously decreasing while p2{x) and p(x) are continuously in- 
creasing functions and p(a) = — pi(a) < 0. We distinguish two cases: 

(1) If p(j) > 0, let 6 G (a, 7] be the only zero of the function p(x), therefore 
Pi{$) = P2(S). The curve Ki+i is obtained from Ki where the parts of Ki 
that determine pi(8) and p 2 {S) are replaced by the vertical line. 

(2) If p(j) < 0, then ^1(7) > ^2(7) and by continuity and monotonicity of the 
function pi(x) there exists only one 5 G (7, 0) which satisfies pi(S) = ^2(7) 
(note that pi((3) = 0). The point 5 together with the areas Pi(8) and ^2(7) 
is marked on the Figure [5j The new curve K i+ i is obtained from Ki by 
replacing those parts of Ki that determine Pi(8) and ^2(7) by a vertical 
line. 



In Figure [6] there is an example of the resulting steps of the above algorithm applied 
to the function given in (25) in time t = 4.25. 




Figure 6. Curves Ki 



0, 1, 2, 3, obtained in three consecutive 



steps of Algorithm 4.1 resulting in the solution at time to. 



We shall briefly describe the method we use to compute the areas needed on 

Let D be a polygon, determined by the points 
(we orient them in counterclockwise direction) and let 



each step of Algorithm |4.1 



T 

Tq(xq, yo) be any point in the plane. Then the signed area of the triangle T TiT i+ i 
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can be computed by 



(23) p ,M+i : = g 



Xi -x yi- y 
x i+ i - x yi+i - yo 

By Green's Theorem one can easily see that the area of the polygon D then equals 



(24) P = $>o,m+i 



i=l 



where T n+ i = T\. 



5. Numerical results 



We have programmed our equal-area method in Mathematica [13] and first com- 
pared the results with the basic Godunov method (which we have also implemented 
in Mathematica) . Further we have compared our method to an advanced Godunov 
method, which is a basis of the widely-used software package Clawpack [2]. Finally, 
we have made a comparison to the very recent software package Particleclaw [9], 
which uses a Langrangean particle method and some information on the charac- 
teristics. Both software packages are freely available on the web. 

The results of these tests are very good. Our algorithm performs favorably both 
in terms of time efficiency and accuracy. The solutions agreed and ours turned out 
to be even more accurate. Figure [7] contains graphs of the solution obtained by 
our method and both above mentioned software packages for the initial condition 



(25) h(x) 



0.9e-* +0.7e-( x ~V + 0.85e"( x+2 ) , xE [-10,10], 
0, otherwise. 




Figure 7. The solution obtained by the equal-area method (thin 
line) compared to Clawpack and Particleclaw, respectively (thick 
dots) . 
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Time complexity of our algorithm depends mostly on finding zeros of a function, 
obtained by the computation of areas of polygons. We used the secant method 
and typically 7 to 12 iterations (9 on average) were needed for 10~ 14 accuracy. The 



number of necessary steps of the Algorithm 4.1 is bounded above by the number 
of stationary points of the function Xt (£) defined in (20). 

We can approximate the error of the position of the shock. Assuming that the 
original curve 7t lies in an e-neighborhood of polygonal approximation line, we 
have an approximation of the area between the "S'-curves" as in Figure [5] 

I e = Ax s 

where I is the length of the "S'-curve", Ax is the displacement of the true shock 
and s is the height of the shock. This gives an approximation of the displacement 
Ax: 

A ■ I 

Ax = £-. 

s 

In Figure [6] polygonal approximation with 1000 points has e less than 6.10 -4 and 
Ax is approximated with 3.10 -3 . The method is quadratical, i.e. doubling the 
number of points would result in decreasing the value of e by factor 4. 

We see the following advantages of our equal-area method. 

• The method is by its nature exactly conservative. 

• The solution is computed for any given fixed time. Hence, the errors do 
not accumulate in time. 

• The method is accurate - the quality of the approximation relies only on 
the quality of the starting approximation of the curve 7j with a polygonal 
line Kq. 

• There is no need for separate treatment of the rarefaction waves, they are 
created on the way where appropriate. 

• The obtained shocks are sharp and propagate with correct speed. Their po- 
sition is obtained automatically by equalizing the appropriate areas. More- 
over, the shock paths are obtained easily by computing the solution for 
some selected times and then simply projecting the shocks from the sur- 
face to the (x, t)-plane (see Figure [8]). 

Acknowledgment. The authors would like to express their gratitude to Ales Zaloznik 
and Marijan Zura for motivation and many helpful discussions. 
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Figure 8. Shock paths in (x,t)-plane for the initial function given 
in p5b. 
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